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We formulate a complete microscopic theory of a coupled pair of bound magnetic polarons, the 
bound-magnetic-polaron molecule (BMPM) in a diluted magnetic semiconductor (DMS) by taking 
into account both a proper two-body nature of the impurity-electron wave function and within the 
general spin-rotation-invariant approach to the electronic states. We also take into account both 
the Heisenberg and the antiferromagnetic kinetic-exchange interactions, as well as the ferromagnetic 
coupling within the common spin BMPM cloud. The thermodynamic fluctuations of the spin cloud 
within the polaron effective Bohr radius of each polaron are taken as Gaussian. 

PACS numbers: 75.10.-b, 75.45.+j, 75.30.Et 
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I. INTRODUCTION 

The bound magnetic polarons, single and molecule, 
represent the charge carrier quantum states formed when 
the magnetization fluctuations influence both their bind- 
ing and a nontrivial thermal behavior. These states have 
been discussed in novel materials encompassing diluted 
magnetic semiconductors,^ ferromagnetic perovskites,- 
and dilute ferromagnetic oxides*^ For example, the origin 
of ferromagnetism in diluted magnetic semiconductors 
(DMS) such as Gai-^Mn^jAs, is intensively discussed in 
the recent years"'^ in view of their potential application in 
spintronics. The origin of ferromagnetismi poses a very 
nontrivial question in view of the dominant role of an- 
tiferromagnetic superexchange in all of those materials 
when the carrier concentration is very lowii Therefore, it 
is crucial to describe these interactions accurately in the 
physically tractable situation. Here we propose a solu- 
ble model of two interacting impurity electrons forming 
a bound magnetic-polaron molecule (BMPM). 

Magnetic interaction between a single electron located 
on a shallow impurity and the localized magnetic mo- 
ments in diluted magnetic semiconductors (DMS) has 
been studied intensively for some time4~— In the first 
period, the influence of classical fluctuations of magneti- 
zation on quantum states of the impurity electron have 
been analyzed^"— , thus leading to the concept of bound 
magnetic polaron (BMP). The main result obtained was 
to demonstrate that thermodynamic fluctuations suffice 
to produce a spontaneous spin splitting of the impurity 
states. A renewed interest in the present decade was 
stimulated by the possibility of ferromagnetic interpo- 
laron interaction^"— which would contribute in a fun- 
damental manner to the origin and properties of ferro- 
magnetic DMS systems and other materialsi^"— In Ref. 
[is ! the authors extended the theory of single BMP— 
to the bipolaron case in the limit when the mutual 
interaction is represented by hopping of electrons be- 
tween the two spatially separated impurities, that leads 
among others to the antiferromagnetic kinetic exchange 
between themJ^ The magnetic cloud is represented then 
by an effective magnetic field, which is oriented arbitrar- 



ily and optimized, but its amplitude is a parameter of 
the approach. In the later version of the approachjii 
the authors approximate the two-impurity wave function 
by piecewise constant values. They consider also a de- 
tailed thermodynamics of the resultant spin model and 
optimize the coupling parameters. The explicit solution 
is discussed in the limit of large interpolaron distance. 
Due to complexity of the interacting BMP pair prob- 
lem existing theories are still incomplete. For instance, 
the model developed by Angelescu and Bhatt^i^ of the 
system of two nonoverlapping polarons is analyzed via 
a generalized Hubbard type Hamiltonian with hopping 
(matrix element t) and Coulomb interaction (energy U) 
turned on, is thought of as being complementary to the 
model of overlapping polaron pair of Ref. ITsl Other 
important and still awaiting for solution problems are 
concerned with the acceptor type interacting BMP's in 
II- VI and III-V DMS. Recent achievements in under- 
standing the ferromagnetism in these classes of materials 
clearly show importance of an accurate treating of the 
non-hydrogeniclike character of the hole wave function 
and its influence by the central cell corrections 

In this paper we introduce and solve BMP molecule 
Hamiltonian for the case of two donors. It is shown that 
our model unifies existing approaches to an interacting 
polaron pair: The Wolff-Bhatt-Durst (thereafter referred 
as to WBD) modelii of overlapping polaron pair and 
the multiple-level generalized Hubbard model of Ange- 
lescu and Bhatt (thereafter referred as to AB) with ran- 
dom fieldsi^. The WBD model is completed by providing 
whole microscopic justification of their Hamiltonian pa- 
rameters. While, the AB model of nonoverlapping, equal- 
magnitude and large polarons is generalized to the case 
of overlapping and arbitrary-magnitude interacting po- 
laron pair, the extension is limited here only to the regu- 
lar Hubbard model (i.e. within one lowest-energy level on 
each polaron site only) . Furthermore, we solve the result- 
ing Hamiltonian within the continuum-medium and the 
effective-mass approximations for the donor case. In this 
manner, we extend the BMP model to the microscopic 
model of BMP molecule. This constitutes an accurate 
ground for the future analysis concerning acceptor-type 
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BMP molecules and the polaron lattices. 

The structure of the paper is as follows. In Section 
II we formulate our model of BMPM. In Section III we 
diagonalize the effective Hamiltonian for the molecule. 
In Section IV the thermodynamical properties are con- 
sidered, whereas Section V contains concluding remarks. 
The Appendixes A and B provide details of analytic cal- 
culations. 



II. THE BMP MOLECULE MODEL 

We start by considering two interacting BMPs in DMS 
with a random distribution of localized spins. Within 
the continuum-medium and the effective-mass approxi- 
mations, one can write Hamiltonian in the form: 



H = Hi + Hn = -Jc S(ri)-Si + S(r2)-S2 
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states. Such an approach greatly simplifies the problem, 
as in practice we have solutions for Hu within the eigen- 
states of {S(Q(, sf^^}. Here we consider more general than 
Heitler-London solution of the hydrogen molecule by in- 
cluding also ionic configurations. In this case, Hu is not 
diagonal in the basis formed from eigenstates of s^^^^ and 
S(Qj. Fortunately, for the off-diagonal elements, the hop- 
ping couples only the different singlet states. 

Next, we express Hj in terms of the creation and anni- 
hilation operators. Details of calculations of the matrix 
representation of Hj within the first quantization for- 
malism are presented in Appendix A. Our methodology 
bases on an explicit demonstration that the postulated 
occupation-number representation of the polaronic part 
of the BMPM Hamiltonian leads to the same expecta- 
tion values as those obtained within the first quantization 
scheme. Thus, we assume Hj in the following form: 
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where m* is the impurity-electron effective mass, a and 
b label the two-impurity sites, I = 1,2 label the two elec- 
trons, Tai, ra2, ni, rf,2 , ri2, ri, r2 and Rat are corre- 
sponding relative distances appearing in the problem, e 
is the static dielectric constant, Jc is the exchange in- 
tegral of the contact Fermi (s-d) interaction between 

localized spins and those of impurity carriers, {s^;}. 

= Si Si5(ri — Ri) is the spin-density operator with 
the sum running over sites occupied by magnetic ions 
(Mn^+). Hjj is the hydrogen-molecule part which sup- 
plemented with the s-d coupling of electrons to localized 
spins {Hi). 

The solution of BMP-pair problem is complexJ^"— 
From one side, a systematic approach to the BMP lat- 
tices and hydrogen molecules can be achieved within 
the second-quantization formalism;^ but from the other, 
the spin part of H is dependent on the positions of the 
localized-moment densites. This fact makes the problem 
more complex to solve. However, it becomes tractable 
when one tries to diagonalize H as consisting of the hy- 
drogeniclike molecule Hamiltonian Hjj influenced by the 
perturbation Hj. This suggests to approach the solution 
of BMPM by selecting as trial states the eigenstates of 
the set of mutually compatible observables {s^g^,s^g^}, 
where Stot is the total spin of the BMPM. In fact, as 
mobile carriers interact via an isotropic Coulomb repul- 
sion, only rotations of the whole molecule conserve their 
relative distance, thus leaving their energy unchanged. 
Moreover, we know that in the Heitler-London approx- 
imation the Hamiltonian commutes with both s^^^ and 
sfot , it is thus diagonal with three- fold degenerate triplet 



where aj^ (flcs) is the creation (annihilation) operator for 
the state on impurity c = a,b with the spin s (s =t, i) 
and the polaron exchange fields Ac are defined as: 
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a, b, (3) 



where M(r) = —g^-B^o^\Xs ^('') Xs/ local 
magnetization per unit volume; A^o = no/v^ is number 
of atoms per unit volume containing fraction x of mag- 
netic atoms, Wa(r) and Wb{v) are the orthogonal molec- 
ular wave functions: 



Wa,b{r) = fi[i^a.b{r) - li^b.air)], 



(4) 



which are build from the single-particle atomic wave func- 
tions and V'fa, being the solution of the corresponding 
single-particle hydrogeniclike Schrodinger equation with 
the effective Bohr radius a^. /3 and 7 are the mixing co- 
efficients: 
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(5) 

determined by imposing the ortogonality and the nor- 
malization conditions, whereas S — (^a(^)IV'b(^)) is the 
overlap integral. It is important to note here, that when 
calculating Ac we replaced the spin density operator S(r) 
by its quantum mechanical average, thereby introducing 
the adiabatic and mean field approximations. This how- 
ever is acceptable here, because the same approximations 
are assumed in the single BMP theory and provide good 
results.^ The BMP's exchange fields Aa and Af, may 
be oriented in an arbitrary direction and we do not make 
any a priori restriction for their magnitude. In our model 
those polaronic exchange fields are overlapping, through 
the presence of the mixing coefficients in the definition 
of Aa and Af,. Next, we introduce the six eigenstates of 
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the carriers total spin operator s = Si + S2: three triplet 
(s = 1) and three singlet (s = 0) states;^ The triplet 
states are: 
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|1> = ^(«lt"H + "Wt) 



|2)-«lXtlO)' |3) = «1XJ0)' 
and the corresponding three singlet states are: 
1 



and 



|4) = ^(«-<«it) 



|5) = ^(«lX; + «IXt)|0)' 



|6) = ^(«IA-4Xt)Io) 
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It is also more convenient to transform the BMP ex- 
change fields Aa and A{, to A+ and according to 
following prescriptions: 

A+ = i(A, + Ab) and A- = i(Aa- Afc), (11) 
which enable us to write the effective polaron fields as: 

A±[M] = — / d^r[\wa{r)\^ ± \wb{r f]Mir). (12) 
This in turn leads to Hj in the form: 

^ ^II(«L(A+ •cr),,,aa«' +aL(A+ •*T),,.a'f„, 

ss' 

+ «L(^" ■ 'T)ss'aas' - al{A- ■ (T)ss'a[,,). (13) 

Now, we select the direction of the global quantization 
axis as aligned with A~. This choice leads to the follow- 
ing 6x6 singlet-triplet matrix representation of Hj : 





-A+COS6I 



A+sine A+C0S6I 
A+s^e'f 

A- 

A- 

\ oA-oy 

(14) 

where 9 and f are respectively, the azimuthal and the 
polar angles between the exchange field A+ and A~. 

For the solution of the BMPM problem, the pola- 
ronic part Hj must be completed with the hydrogenlike 



molecule part Hu, written also in the second quantiza- 
tion form. Such representation of Hu is knowni&^~ so, 
we can write it directly as: 

Hii = eaUa + ebnb+t'^ (aL^as + aL^fc^) +Uana^nai 



+ Ubm^Ubi - 2 Jsi • S2 -I- yK + j "a^f" 
-t- J {al^al^abiUb^ + H.c.) + V^[{nas +"-„«) 

S 

X {al-abs + al-aas)], (15) 

where, Cc is the atomic level position, t the hopping 
between the impurity states and Uc the intra-atomic 
Coulomb interaction. The remaining terms represent re- 
spectively, the interatomic (Heisenberg) exchange, the 
pair hopping, and the so-called correlated hopping. In 
the singlet-triplet basis, Hu has relatively simple form^^, 
so we may write: 
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(16) 



where: 

L = 2e + K~j, Le = 2e + U~j, andT = 2{t + V). (17) 

For identical impurities both the atomic level positions 
and Cb and the intra-atomic Coulomb interactions Ua and 
Ub are equal (i.e. = = e and Ua = Ub — U). We note 
here, that Hu is solved in the local coordinate system, 
with z axis aligned with the direction of the molecular 
bond, and then Hamiltonian Hu is reoriented towards 
the global quantization axis. This last step is quite trivial 
due to the rotational invariance of Hu. 

Having determined the BMP molecule Hamiltonian, 
we discuss next its relation to the with existing mod- 
els of Angelescu-Bhatt (AB) and Wolff-Bhatt-Durst 
(WBD)J^~— These earlier models of interacting BMP 
pair, while being significant theoretical achievements, are 
subject to certain important limitations, clearly stated 
by their authors. Here, for a purpose of completeness, 
we enumerate some of them. As in our model each site 
contributes a single s-type orbital, the AB model ap- 
pears as more advanced, because it accounts for excited- 
states transitions from Is level up to 3d level. Never- 
theless, it should be kept in mind, that AB model treats 
BMPs as nonoverlapping, equal and large field within 
the perturbatoion approach with respect to their mu- 
tual interaction, hence neglecting important processes 
even in the simplest one-level situation. For instance, 
AB model does not predict ferromagnetic ground state 
in that case. Their assumption about equal magnitudes 
of BMPs exchange fields creates a serious problem with 
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FIG. 1. (Color online) Mixing coefficients 7 and 7^ as a func- 
tion of interimpurity distance Rab expressed in the units of 
as- 




RabEae] 



FIG. 2. (Color online) K, and y^K (in units of a) as a 
function of interimpurity distance Rat expressed in the units 
of a_B. 



treating them as thermodynamically fluctuating quan- 
tities. From the other side, trusting in predictions from 
WBD model is weakened by not fully microscopic deriva- 
tion of their Hamiltionian, including their parameters. 

Therefore, it is important to recognize consequences of 
such approximations, even in a relatively simplest case. 
In this work the model is formulated as applicable beyond 
those limits namely, our BMPs are overlapping, have ar- 
bitrary magnitude of the exchange fields and their mutual 
interaction is accounted for within molecular formalism. 
Moreover, we are able to consider also thermodynamics, 
including thermal fluctuations of BMPs exchange fields. 

To demonstrate the generality of our approach more 



clearly, we map polaronic part Hj into corresponding 
parts of the WBD model Hamiltonian (the whole WBD 
model contains also the antiferromagnetic Heisenberg ex- 
change interaction between carrier spins, which we do not 
specify explicitly). Namely, we can rewrite Hj in the spin 
operator form: 



Hr 



Aa • Si + Af, -82, 



(18) 



and take into account the explicit form of the exchange 
fields Aa and which yields: 



9^J'B 
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This expression should be compared with following pola- 
ronic part of the WBD model Hamiltonian)ii 

HwBD - K[si • Si + ?2 • S2] + K'{si + S2) ■ S3, (20) 

where now K and K' are the WBD model parameters and 
S3 is the total spin of magnetic ions in overlapping region. 
A direct mapping can be established by introducing the 
following definitions: 



K = aP\ K' = —fK, Si(2) 



and S3 = / i^*a^Pb + ^b^a) 
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(21) 



which allow to rewrite Hj into the final form: 

Hj = K[si • Si + S2 • S2] + K'{si + S2) • S3 

-|-72i^[si-S2+S2-Si]. (22) 

In Fig. 1 we have plotted 7 and 7^ as a function of inter- 
polaron distance Rat, where in Fig. 2 we show the same 
radial dependence for K, jK and j^K. As can be seen 
from these plots only K is different from zero at large dis- 
tances and in this limit is equal to that corresponding to 
isolated BMPs. In Hj, the parameters 7 and 7^ appear as 
interaction couplings and the corresponding three terms 
describe contribution from the zero-, first- and second- 
order processes, respectively. One can also see, that at 
large interpolaron distances, for which 7^ — > 0, Hj can be 
reduced to the bipolaron part of WBD model Hamilto- 
nian. In this manner, we have provided the microscopic 
derivation of the WBD model and the microscopic mean- 
ing of their parameters. Also, the derived relation of 
BMP exchange fields A^, and Ah to the carrier ortognal- 
ized single particle wave functions shows that the mag- 
nitude of these quantities depends on the interpolaron 
distance Rab, the feature which is completely neglected 
in the AB model. 
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III. SOLUTION OF THE BMP MOLECULE 
HAMILTONIAN 

Diagonalization of the BMPM Hamiltonian requires 
the sohition of the eigenequation HV = EV, with H 

I 



H = 



Its form does not allow for an exact analytical diagonal- 
ization. Therefore, before solving it numerically, we dis- 
cuss first some limiting situations, for which approximate 
analytical solutions can be obtained. 

A. Saturation limit 

We first solve BMPM effective Hamiltonian for pair of 
donors in the saturation limit, i.e. with M(r) = Msat- 
In this case, the solution is simplified greatly within our 
approach, because we can write: 

a 

A+ = Msat = ?>/2aN and A_ = 0, (24) 

where N is the number of Mn ions within BMPM which 
contribute to magnetization. Subsequent diagonalization 
of H leads to following exact solution for the enigenvalues 
Ei, i = 1,...,6: 

/ ^, \ 

L + 5/2aN 

L - 5/2aN 

i [(i + Le) + 4j + ,J{L-L,f+AT^] ■ (25) 

+ Le) + 4j - ^J {L - Lef + AT^] 
\ Le J 

The eigenvalues E^^^ contain antiferromagnetic kinetic 
exchange interactional in an explicit form, that com- 
petes with the direct Heisenberg exchange. Nature of 
the ground state depend on the sign of the expression 
AE = £'5(A_) — £'3(A+), which determines the domi- 
nant exchange interaction aligning individual polaron po- 
larization clouds;for positive value the ground state is fer- 
romagnetic. Note that states belonging to the eigenvalus 
El, E2 and E3 are the triplet states. We plot in Fig. 
3 the singlet-triplet splitting AE as a function of inter- 
impurity distance Rab, calculated for the material pa- 
rameters corresponding to Cdo.95Mno.01Se; Nqu = 0.28 
eV (line a)M' This naturally is only a formal dependence, 



given by: 



(23) 



r 




FIG. 3. (Color online) Exchange splitting AE dependence 
on inter-impurity distance Rab for BMPM with effective con- 
centration X = 0.027 of magnetic ions. Lines a and 6 
are calculated with material parameters corresponding to 
Cdi-^Mn^Se and Gai_a;Mna;As, respectively, within the s- 
state approximation for the impurity single-particle wave 
function. 



which shows a potential strength of a force stabilizing 
ferromagnetic state. One sees also that if the two impu- 
rities are to close each other, then nonmagnetic hydrogen- 
molecule spin-singlet ground state becomes dominant. In 
s-type II- VI DMS, the magnetic susceptibility x of local- 
ized magnetic ions at low temperatures takes the form 
of the Curie- Weiss law Cm/{T + Tq), with To > 0-^ 
Therefore, the magnetization is certainly not saturated 
and considered case appears as describing a non-realistic 
situation. The situation changes for p-type ferromag- 
netic III-V DMS in which magnetization saturates and 
presented here model provides a quantitative argument 
that for certain range of interpolaron distances the in- 
teractions between BMPs can stabilize ferromagnetism. 



L ^-Hs^giv ^+s^g-»v A- 

A+si^e-'^P L-f A+C0S6' 

A+SinO^^,p L-A+COS6I 

A- L + 2j T 

T Le + 2j A- 
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In Fig. 3. the line (b) was calculated with values of 
parameters corresponding to p-d exchange Nq/S = —1.2 
eV confirming our last sentence. Nevertheless, it should 
be stressed that the nature of holes wave function in 
zinc-blende semiconductors is much more complicated, 
then the used here simple s-typc donor wave function. 
In other words, the acceptor BMPM is not quite correct 
theory.' ~' We postpone this discussion to a separate 
publication. 

B. Lcirge inter-polaron distance limit 

It can be derived by making the transformation 
R~^HR on the whole 6x6 Hamiltonian matrix, with 



R being the rotational matrix which diagonalizcs Hu. 
In explicit form R is the sparse matrix, with all el- 
ements equal zero, except for R44 = i?55 = cos(i5), 
R45 = — sin((5) and R54 = sin(^), and the rotation angle 
5 is defined through: 



ta„(2*) = ip:n. 



This leads to H in the form: 
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2sin(e) 



2s^gi^^+ 2si^g_i^^+ ^-cos(5) -A-sin(^) 
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(27) 



where the diagonal elements L4 and L5 are: 

i4,5 = 2ea + l/2{K+U)+j±l/2[{U-Kf+lQ{t+Vff'^. 

(28) 

In the considered here limit of large interpolaron dis- 
tances one can neglect the terms proportional to sin((5). 
In practice, this limits us to the region with Rah Jc^ 3.5 
as. In that limit, we can write the resulting sixth-order 
eigenequation as: 

{E^ - E{L^ + Lq) + L^Lq - A-'^) 

{[{E - Lf - A+2](i;2 - E{L + U) + LL4, - A-2) 

- (A+ X A-)^} = 0. (29) 

The eigenvalues E^ and Eq can be easily determined, but 
presence of the vector product (A"*" x A~)^ still compli- 
cates calculations of remaining four eigenvalues in simple 
terms. Formally, this equation can be solved exactly, as 
in a principle all roots of the foiir order equation can 
be found analytically. However, solutions will have com- 
plicated analytical form, which precludes further explicit 
analysis. To overcome this difficulty we regard the ex- 
pression in {...} as a function of £, say F{E), and write: 



where 



F{E,e) = G{E) + A{0), 



A{e) = -{A+y{A-y sin'ie), 



(30) 



(31) 



the definition of the function G{E) is self-explanatory. 
Now, it is easy to observe that the quantity A{0), which is 



independent of the cncirgy E. influences F{E, 9) through 
a downward shift of G{E). Thus we need to analyze prop- 
erly G{E) and at least take into account approximately 
the presence of A{9). Therefore, we can expand F{E,6) 
in terms of Taylor series for each eigenvalue separately 
around the zeros of G{E): 



F{E,0)=A{e) + 



1 dG{E) 
T\ dE I, 



{E — Eoi 



+ 



1 d^G{E) 



{E - Eoif + ... = 0. (32) 



2! dE^ u 

We solve those equations to the first order and obtain: 
Ei = Eoi- ^ sG^^O ' for i = 1,...,4, (33) 

TT dE 

In effect, we find the following six zero-order eigenvalues 



( 1/2{L + Li~ v'4(A-)2 + (L - Lif) \ 
L + A+ 

L-A+ 

1/2{L + L4 + v/ 4(A- )2 + (L-L4)2) 
^5- V4(A^ 



1/2(^6 



(34) 



V 1/2(^6 + L5 + V4(A-)2 + (L6-i5)2) / 



Analysis of this expression shows that the nature of the 
ground state is determined by the sign of AE = E^ — E3 
expressing the difference between parallel and antipar- 
allel configuration of the exchange fields. Moreover, for 
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A~ = 0, the eigevalues have the same functional form 
as in the exact solution, except that now A+is not cal- 
culated for saturated situation. 



1. Asymptotic solution for Rat — > oo 

In the limit Rab — > oo, we can neglect the terms de- 
scribing interactions between BMPs and two ionic con- 
figurations. Naturally, in resulting expression and 
cannot be neglected. Next, we diagonalize Hji^^ 



i.e. write: 

Hr^,^^^V^EV. (35) 

This equation leads to the related fourth-order equation 
of the form: 

[i2E^ - Ef - (A+)2][((2i?o, - Ef - (A-)2] 

- (A+)2(A-)2sin2(0) = 0, (36) 

where Eoo — — ^b- The corresponding four eigenvalues 
are: 
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(37) 



Our asymptotic solutions still depend formally on the 
angle d between A+and A~, whereas the corresponding 
solutions for the two isolated BMPs: 

/2£;o,-|(A„-)-Ab)\ 

- 2(Aa - Ab) , , 

2i5oo + |(A,-Ab) ^^"^ 
V2£;oo-(-|(A„ + Ab) / 

are clearly free of such angular dependence, what reflects 
an uncorrelated character of the spatial orientation of the 
polaron fields. In Appendix B we prove the equivalence 
of our solutions given by Eq. (37), with that for the 
two isolated polarons. Note, that for finite interpolaron 
distances, the limit A — > leads also to a correct solution 
of the hydrogeniclike molecule. 

IV. THERMODYNAMICS 

Thermodynamic fluctuations of magnetization may 
strongly influence behavior of the system. Having de- 
termined the eigenvalues, we can construct the free en- 
ergy of the BMP pair. To determine these properties 
we extend the previous approach^ devised for a single 
BMP. For the case of single BMP, the thermodynam- 
ics has been derived by including the contribution com- 
ing from localized magnetic moments starting form the 
Ginzburg-Landau Hamiltonian: 




where the two phenomenological parameters, k and x~^j 
are the exchange stiffness constant and the inverse static 
susceptibility, respectively. For this case, the probability 
distribution of the exchange field has been found in the 
form: 

P(A) = y DM(r)P[M(r)](5(A - A[M(r)]), (40) 

where the probability distribution of M(r) is defined as 

P[M(r)] = Ce^'^^^e'^^^ , (41) 

and C is a normalization constant. The functional in- 
tegration in Eq. PO)) expresses a summation of the con- 
tributions coming from all space profiles of magnetiza- 
tion {M(r)} contributing to given value of the exchange 
field A. This necessitates the functional integration 
over all possible "paths" of {M(r)}, with the proba- 
bility density P[M].^ Generalization of Eqs.(l39l)-(|40l) to 
iV -polaron case leads to the following A^— component 
Ginzburg-Landau Hamiltonian: 



i/s[Mi,...,Mjv] = 




Next, we define the probability distribution of exchange 
fields Ai, with i — 1,...,N, in the following form 



P(Ai,...,Aw) = 

N 

C / l[DM,e'^"''^ S{A, - A[M„u;,])e-'^^^- (43) 



8 



A+, A- is: 



0.20 



0.15 



0.10 



0.05 




FIG. 4. (Color online) Interpolaron distance dependence of 
the functions f+{Rab) and f-{Rab) describing the fluctuating 
exchange field, as defined in Eq. (I47|) . 



where AFt is the electronic part of the total free energy 
of N polarons system, C" is now the overall normalization 
factor, and 



M,(r)|?«,(r)|'dV (44) 



Transformation made in Eq. (1431) defines our dynamic 
variables and the functional integrations in this equa- 
tion can be carried out for each A.[M.^,Wi] with the 
methodology developed in Ref. 8. 

When coupling between BMPs is important, our eigen- 
values are given by Eq. ([55]) and the two-isolated BMP 
problem transforms now to the problem of interacting 
polarons with A+ and A^. 

Thermodynamics of this system can be derived from 
Eq. which after executing the functional integra- 

tions provides the physical free energy: 



F = kBT\n // d'^A+d^A-p(A+, A"), 



(45) 




knT . 



, (46) 



and the eigenvalues Ei arc the roots of Eq. p9| , whereas 
the two new parameters, e+ and are defined as fol- 
lows: 



4 (ffMis)' 



d\[\wa{r)\^ ±\wb{r)\ 



2l2j3 



(47) 



where the probability distribution of the exchange fields 



Note that P(A^, A+) has a correct asymptotic behavior 
for the interpolaron distance Rab — >■ oo. We stress also, 
that the parameters e+ and £_ are material dependent 
{x is the spin system susceptibility) and correspond to 
the single-polaron parameter £p of DS.^ Their Rab de- 
pendence is crucial, as it affects the magnitudes of either 
A+ or A^. In Fig. 4 we plot this dependence. One 
sees, that for all finite distances /+ > /_, i.e. the field 
responsible for the triplet configurations of the two impu- 
rity electrons is dominant. Note also, that even though 
the parameter /3 — >■ c» in the limit Rab 0, the corre- 
sponding functions f+{Rab) and f-{Rab) are finite then 
(see also below). 

The dominating character of £+ may result in the 
ferromagnetic-ground-state appearance of the BMPM. 
However, to prove that explicitly we have to include 
also the effect of the exchange fields on the eigenvalues 
{Ei}. This can be done by determining the most prob- 
able values of the fields A"'' and A~. For that purpose, 
we consider first the lowest-order solution assuming that 
Ei « Eoi. In analogy to the single-polaron theory, we 
can now calculate the most probable values A+ and A~ 
from the conditions: 

[9P(A+, A-)/dA+]^ = [aP(A+, A-)/9A-]^ = 0, 

(48) 

which lead to the following system of two coupled equa- 
tions: 



(1 + A + B)A- -4e_ 



^\j' + {A-y 



zC tanh 



IJ^ + (A-)^ 



8(1 + A + B)A-e-kBT = 0, (49) 



and 



{B-^ + 1 + AB'^)A+^ - 4e+ A+^ tanh 



knT 



8(1 + B-^ + AB-^)A+e+kBT = 0, 



(50) 
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where A = A{T), B = B{T) and C(T) are defined as: 



A{T) = cosh 



exp 



J6 
2 



2 



cosh 



-1 



(A- 



(51) 



B{T) = cosh 



" A+ " 






exp 



L 



cosh ^ 



iJ^ + (A-)^ 
knT 



(52) 



and 



C(T) = 1 + 





V'iJ| + (A-)^ 


v/iJ^ + (A-)^^''^^ 


feBT 


^i-^i + (^)^ sinh 







(53) 



with J = Lq — L5 and J4 = L4 — L. This system of 
coupled equations leads to the following three pairs of 
the solutions at T = 0, depending on the sign of AE = 
E5 — E3: 

(i) A^' = {2e-f - {J/2f, and A+ = 2e+, 

for AE = 0, 



A+ = 4e+, and A" 



0, 



for AE > 0, (54) 



and 



{in) A+ = and A" = V(4eZy 



- (^/2)^ 

for AE < 0. 



The condition AE{Rab) = defines the critical distance 
Rc for Rab > Rc the ground state is ferromagnetic. 
Therefore, the driving interaction which aligns at T = 
individual polaron polarization clouds can be defined at 
Rc in the form AEc = E^{iii) — E^lii). Let us then dis- 
cuss explicitly the limit T ^ 0, for which the ground 
state is determined by the sign of the expression: 



•^(^«'') . L,-J\j{R,bf + {A-{R^b)y 



L-A+{Rab) =f{Rab). (55) 



This expression for negative values yields the ground 
state belonging to eigenvalue E^, i.e. a mixture of two 
singlet states of impurity electrons, whereas for f{Rab) > 
it leads to ferromagnetic, (sj^j = ±1) ground state with 
the corresponding eigenvalue E^. Therefore, the condi- 
tion f{Rab) = defines a critical interpolaron distance 



Rc, at which the crossover from magnetic to nonmag- 
netic ground states occurs. A direct analysis shows that 
limT-s-o A{T) = and limx->o C'(T) = 1, whereas an im- 
portant role on the cliaractcr of the ground state is played 
by the function B{T). This become clearly visible if one 
determines its value for T ^ 0. Namely, B{T — > 0) ^ 00 
for f{Rab) > 0, so S-i(T ^ 0) 0. It can be readily 
seen, that in this case we have: 



A+ = 4£. 



and A" 



0. 



In the other case, i.e. for f{Rab) < 0, B{T 0) 
and we find: 



A- = V(4£-)' - {J{Rab)/'^?, and A+ = 0. 



(56) 


(57) 

Such peculiar behavior of B{T — >■ 0) stabilizes the ground 
state and moreover, makes our analytical solution exact, 
because Ei = Eoi under these circumstances. For com- 
pleteness, we need to determine still B{T — >■ 0) at Rc- It 
can be shown that B(T 0) 1, with: 



A- = V(2£-)2- (J(i?c)/2)2, and A+ = 2e+. (58) 

Even at low nonzero temperature solutions for A+ and 

A^ become complicated and temperature dependent. 
Therefore, considering behavior of the state population 
probabilities at low T, one can propose an approximate 
expression for determining Rc'. 



/+ - bf- 



L-6L5 



hJ 
4^ 



(4x/_)2 



w{Rab), (59) 



where Jn-iiRab) = J{Rab) is taken for the hydrogen 
molecule (i.e. e = 1 and m* /rric = l,me being the elec- 
tron mass), whereas Rc is expressed in units of ob, and 
K is defined as: 



a Y e me 

H = 

m* 



(60) 



and contains all the key material parameters that deter- 
mine Rc, provided that the parameters L, L5, J15 and J56 
are calculated for the H2 molecule. Note that neglecting 
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FIG. 5. (Color online) Plot of w{Rab), defined in Eq. (59), 
for three DMS characterized hy >c — fcxo, where the value xo 
corresponds to that for Cdo.95Mno.05Se and fc > 1. 



FIG. 6. (Color online) Temperature dependence of average 
exchange fields A"*' and A~ for DMS characterized by = 
25x0. 



the quadratic term in the condition w{Rab) = and for 
6 = 1 is equivalent to AE = 0. Next, we discuss the Rc 
dependence on material parameters. 

In Fig. 5 we plot w{Rab) with b = 1.1 (T ~ I.IK) 
for DMS characterized hy x = k>to, with k — 4,5, 
and 15, where the value >fo corresponds to that for 
Cdo.95Mno.05Se (e = 9.4, m*/me = 0.13, (A = 2ep < 1 
meV)).— One sees from Fig. 2 that for Cdo.95Mno.05Se, 
w{Rab) < for all Rat', then the Sj^j = spin configura- 
tion of the impurity electrons is always stable. 

So far we have treated the magnetic susceptibility 
X of localized magnetic ions as a material parame- 
ter neglecting its strong temperature dependence. In 
Cdo.95Mno.05Se and at low temperature, x takes the 
forme of Curie- Weiss law x = Cm/{T + To) with To = 1.2 
K. Properties of the system at low temperatures may be 
deduced from analysis of the averages of A+ and A~ , as 
they reflect properties of the most probable values A^ 
and A expressed via conditions (i) — (iii). In order to 
see that within our approach the ferromagnetic config- 
uration of the impurity electron is possible in general, 
we have displayed in Fig. 5 the numerical average of 
A"*" and A~ vs. T for exemplary value of the parameter 
>c{T) = kKQ{T), with k = 25. Such value of k for n-type 
DMS is naturally outside of accessible range of material 
parameter, but for the p-type DMS the p-d exchange is 
about 4—5 times stronger then that for the n-type, then 
such high value of k is justified. For this value of k one 
finds AE'c « 0.47 meV (calculated at 2 K, with Tq = 1.2 
K and for Rab — 3.65 &b), which allow for a very rough 
estimation of the Curie temperature Tq, in the mean field 
approximation for a 3D cubic BMP lattice, would be 22 
K. 

For completeness, we plot in Fig. 6 temperature de- 
pendence of the average values of A+ and A~ taking the 



parameters for Cdo.95Mno.05Se at Rab — 2 and 4.25 a^. 
As can see, the numerical results confirm that for this 
material a nonmagnetic (sfoj = 0) configuration of the 
impurity electrons is always stable. 



V. CONCLUSIONS 

In this work we have developed the model of BMP 
molecule (BMPM) consisting of two overlapping po- 
larons, which have arbitrarily oriented in space their po- 
laronic spin clouds of arbitrary magnitude. Their mu- 
tual interaction is accounted for within the molecular 
electronic states. We also succeeded in deriving ther- 
modynamics of these states that allows to analyse their 
properties at nonzero temperature. These features allow 
for a better understanding of the earlier results, includ- 
ing the range of their applicability. More directly, an 
important role of an accurate inclusion of the impurity 
electron wave functions has been demonstrated. What is 
also important, the microscopic interpretation of Wolff- 
Bhatt-Durst model Hamiltonian-parametersi^ has been 
provided. Simultaneously, the regular case of the many- 
level generalized Hubbard model with random fields of 
Angelescu-Bhatt^^ of nonoverlapping, large polarons has 
been extended to the case of overlapping polarons form- 
ing the BMP molecule. 

Generally, our model of BMPM in DMS confirms a 
possibility of the ferromagnetic ground-state appearance 
for certain DMS materials and appropriate distances be- 
tween polarons. One of the important and new results 
of our approach is a direct incorporation of material 
parameters into the BMP molecule model. Presented 
here numerical calculation for BMP molecule embodied 
in Cdo.95Mno.05Se shows, that for all interpolaron dis- 
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tances the polaronic molecule ground state spin config- 
uration is S(j,( = for the impurity electrons. Then, in 
this DMS a non-magnetic ground state of BMP molecule 
is always stable. From the other side, our model pre- 
dicts a ferromagnetic (spin-triplet) ground state of BMP 
molecule for material parameters corresponding to p-type 
DMSs. The nature of the ground state depends on mag- 
nitude of the parameter >f, defined by Eq. (1601) . which 
is proportional to the magnetic susceptibility of the spin 
system of host DMS. Therefore, the RKKY interaction 
responsible for the low temperature ferromagnetic order- 
ing in p-type DMSs enhances ferromagnetic ground state 
of BMPM when approaching the critical region from the 
high-temperature side. However, the results for the p- 
type DMSs must be analyzed within the approach, which 
accounts correctly for the non-s-type character of BMP 
carrier wave function, as well as non-Gaussian contribu- 
tions to the free energy functional. Such formulation of 
BMPM is under consideration and will be presented sep- 
arately. 



Appendix A: The matrix representation of Hi 
within first quantization formalism 

We assume that the exchange field acting on impurity 
electrons is relatively weak and that orbital moments can 
be disregarded. As a consequence, the two-impurity elec- 
tron wave function may by written as a product of func- 
tions, which depend on the spatial and the spin variables 
separately i.e. it is of the form 

^c(|Si| ;ri,r2,cri,cr2) = Xs Xo.i^2*c(ri, , 

(Al) 

where x |Sj| is the wave function of the spins, and Xaia2 
is the spin part of that for the carriers, whose spatial 
wave function is 5'c(ri,r2), where c runs over all pairs 
{ab,aa,bb). We take the two-electron functions as com- 
posed of orthogonal molecular wave functions Wa,b{i'i)- 



(A2) 



contain Is-hydrogenic-like wave functions 'i/'a ^'^^ 
with the effective Bohr radius a^jwith /3 and 7 being 
the mixing coefficients. These orthogonal single-particle 
wave functions are used to construct the two-particle six- 
dimensional space. The space part is identical with that 
for the hydrogeniclike molecule. Therefore, we can write 
down directly the following six two-electron wave func- 
tions: 



*l,4Xl. 



1 



V2 



Wairi)wb{r2) T Wa(r2)wb(ri)] 



for the spin-triplet state with s^^^ = and the singlet 
state respectively and 

^2.3X2,3 = ■^['Wa{ri)wb{r2) - u;Q(r2)u;6(ri)] 



X [x±i/2(l)x±i/2(2)] , (A4) 
for the remaining triplet states with s^^^ il, and 
1 



1 



X ^[Xi/2(l)X-i/2(2) ± X-i/2(l)Xi/2(2)], (A5) 

for the two ionic singlet states. 

Next, we find the matrix representation of Hi by av- 
eraging over the spatial coordinates, but keeping out the 
carrier spin still in the operator form: 

(^-^ \Hi\ - — [si • / *;M(ri)*„d3ridV2 

+ S2- j ^*JA{Y2)^md\ld\2]. (A6) 

Taking into account the symmetry properties of spatial 
wave function with respect to the transposition of the 
carrier coordinates, the matrix elements of Hi can be 
recast to the form: 

{^n \H\ = Si • Ai[M,«'„, + S2 • A2[M,*„,^'„] 



= (si + sign{^,,)sign{-^„,)s2) ■ A[M,«'„, , (A7) 

where the factor sign{...) = ±1 and expresses the parity 
of the corresponding wave function. The exchange fields, 
A;[M,vI/„, vI/,„] are thus defined as: 

Az[M,*„,^„] EE — / / *;M(r,)*„dVic;V2, 
5Mi3 J J 

(A8) 

Note that the fields Ai[M,^'„, ^Pm] are equal to 
A[M,^'„, except for the sign. An explicit calcu- 
lation provides a remarkable reduction of the number of 
the exchange fields. Namely, only the following two ap- 
pear in the final BMP pair Hamiltonian matrix: 



A±[M] = — [ S 



r[\wa{T)\' ±\wb{r)nM.{T) . (A9) 



X -^[xi/2(l)X-i/2(2) ± X-i/2(l)Xi/2(2)] , (A3) 



The labeling of the exchange fields A+and A^ originates 
from the wave function parities (e.g. corresponds to 
the triplet-triplet matrix elements and "-" to the singlet- 
triplet ones). 

Here we outline the details of calculations of the ma- 
trix elements ^(j)**"*'"*"* \Hbmp\ (1)**°*'*'°'^ in the singlet- 
triplet basis. According to the definition given in the 
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text : = ^'^Xt*°*'*'°' , and taking into account that 

Hamiltonian Hu does not contain explicitly the carrier 
spin variables, one can write: 



Generally, one can assume that in the selected coor- 
dinate system direction of A+[M] is described by the 
angles and (j). Therefore, evaluation of Eq. (AlO) with 
these exchange fields is the simplest in the appropriately 
rotated spin basis, namely: 



E^5rMxr'°''"'°' |(si + sign{Tv)s2) ■ A^"^]] x^*-*''*"* 



(AlO) 



where Si = Si (g) 1, S2 = 1 S2 are the carrier spin opera- 
tors and the corresponding spinors ^r*"*'**"* are explicitly 
defined below: 



1 



XT = 7|[^i/2(l) ® X-l/2(2) - X-l/2(l) ® Xl/2(2)], 

(All) 



1,0 



;^[Xi/2(l) <^ X-i/2(2) + X-i/2(l) <^ Xi/2(2)], 



1,1 



Xl/2(1) ® Xl/2(2), 



(A12) 
(A13) 



and 



X-'"'=X-i/2(l)®X-i/2(2). (A14) 



To short-hand notation of the tensor product symbol are 

omitted below. Wc recall that we have fixed the direc- 
tion of the spin quantization axis as parallel to the local 
exchange field A~[M]. 

I 



Xi/2 = X00e'2 cos - - 



X-i/2 = Xe0e '2 sin- 



e'^ sin 



2' 



e '2 cos-. 



(A15) 



(A16) 



where the angles and (p describe the spin state with an 
arbitrary direction opposite to that defined by and <p. 

Prom Eq. (AlO), the diagonal matrix elements can be 
written as follows 



xP''''''-'\Hbmp\x>*''' 

= Er + (x?" I (si + S2) • A- 1 xr*"* ) ■ (A17) 



Note that {xT\Hbmp\x+° 



Ej^ because x+° 



transform as a scalar under rotations, i.e. 



Xl/2(l)X-l/2(2)-X-l/2(l)Xl/2(2) 

= Xe0(l)X|^(2) - X5^(1)X£42). (A18) 



Also, the matrix element ^X-° \Hbmp \ X-° ) can be com- 
puted in the following manner 



xi'°|(si+S2)-A--|xi'° 

= ^ E (xi/2(l)X-i/2(2) + X-i/2(l)Xi/2(2) \% ■ A— I Xl/2(l)X-l/2(2) + X-l/2(l)Xl/2(2) 



1/2,-1/2 



1/2,-1/2 



= ^E E (x.(l)x.'(2)|s.-A--|x.(l)x.'(2)) = ^l] {x.mi-^--\xM 

1,2 ^24>2^i-2-^'i'2 

= 9 E E (x-(^) A-- = -A-- [(cos2(0 - sin2(0) + (sin2(0 - cos'ii))] = 0, (A19) 



and thus (x-° \Hbmp\ X-°) = E- 
The other matrix elements are: 



X-^ \Hbmp\X-^) = 



E (^1/2 (*) I Si • A+ 1 Xi/2 («) ) = A+ cos ^, ( A20) 



and 



i.-i I tj I i>-i 
X- \Hbmp\X- 



E(^-i/2W|^i- A+IX-1/2W) =-A+cos0. (A21) 



Finally, the non vanishing off-diagonal matrix elements 
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are: 



X+ \Hbmp\X- 



A-(xrix°/). (A22) 



where, Lp and are respectively, the azimuthal and the 
polar angles between the exchange field A+ and A~. 
This form of the Hamiltonian matrix is identical to that 
taken in Sec. II for a detailed analysis. 



Appendix B: Equivalence of Eq. (35) with the 
solution for two isolated BMPs 



\Hbmp\ X-'7 = {X-"" |(si + S2) • A+l X 

= ~7nYl (x-i/2{j) |sj • A+i xi/2{j) 



l,l\ _ / 1,0 



+ 1 



V2- 



= —A+e''^sm9, (A23) 
V2 



and 



X-° \Hbmp\X--') = (x-° \Hbmp\ X-'') • (A24) 



Next, we select a convenient direction of the quantiza- 
tion axis for the BMP-pair problem, which is selected as 
aligned with the field A~. In effect, the matrix repre- 
sentation of Hj in the six dimensional basis of the two- 
electron wave functions is: 



^+sine_^zv ^+sme_^-^.f A" \ 



V2 
V2 

A- 







-A+ cos( 












A- 

A- y 

(A25) 



In this Appendix we demonstrate that our solution of 
the asymptotic case for Rat — > oo is equivalent to that 
know for two isolated BMP, hence we prove the equiva- 
lency between Eq. (37) and Eq. (38). To do this we need 
calculate first: 



2A+A" cos(6i) ^ (A+ + A-f - {A" 
Taking now into account that: 

(A+ + A_)2 = AAl 

and 



(A- 



one can find: 



(A 



+ ^2 



(A-y = 2iAl + Al), 



A+A" cos(6i) 



(Bl) 



(B2) 



(B3) 



(B4) 



Substitution of Eqs.:(B3)-(B4) into Eq.(37) gets finally: 



2E^ 



+ i(A, + Ah) \ 
-|(A„ + Ab) 



(Aa - Ab) 



(B5) 



\2E^ - ^{Aa - Ab) ) 



These eigenvalues can be found also taking independent 
linear combination of the solutions for the two isolated 
polarons. This constitutes the equivalence between the 
two solutions. 
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